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Star patterns, reminiscent of a wide range of diffusively controlled growth forms from snowfiakes 
to Saffman- Taylor fingers, are ubiquitous features of ice covered lakes. Despite the commonality 
and beauty of these "lake stars" the underlying physical processes that produce them have not been 
explained in a coherent theoretical framework. Here we describe a simple mathematical model that 
captures the principal features of lake-star formation; radial fingers of (relatively warm) water-rich 
regions grow from a central source and evolve through a competition between thermal and porous 
media flow effects in a saturated snow layer covering the lake. The number of star arms emerges 
from a stability analysis of this competition and the qualitative features of this meter-scale natural 
phenomena are captured in laboratory experiments. 
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I. INTRODUCTION 

The scientific study of the problems of growth and form occupies an anomalously broad set of disciplines. Whether 
the emergent patterns are physical or biological in origin, their quantitative description presents many challenging 
and compelling issues in, for example, applied mathematics biophysics [2], condensed matter Q and geophysics 
U wherein the motion of free boundaries is of central interest. In all such settings a principal goal is to predict the 
evolution of a boundary that is often under the influence of an instability. Here we study a novel variant of such a 
situation that occurs naturally on the frozen surfaces of lakes. 

Lakes commonly freeze during a snowfall. When a hole forms in the ice cover, relatively warm lake water will flow 
through it and hence through the snow layer. In the process of flowing through and melting the snow this warm water 
creates dark regions. The pattern so produced looks star-like (see Figured]) and we refer to it as a "lake star" . These 
compelling features have been described qualitatively a number of times (e.g. 0, @, 0]) but work on the formation 
process itself has been solely heuristic. Knight [5| outlines a number of the physical ideas relevant to the process, but 
does not translate them into a predictive framework to model field observations. Knight's main idea is that locations 




FIG. 1: Typical lake star patterns. 
Connecticut, 8 March, 2006. 



The branched arms are approximately 1 m in length. Quonnipaug Lake, Guilford, 
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FIG. 2: Schematic of the geometry of the model. The perspective is looking down on a nascent star. The equations (refer to 
text for numbering) are shown in the domains of the system where they are applicable. 



with faster flow rates melt preferentially, leading to even faster flow rates and therefore to an instability that results 
in fingers. This idea has features that resemble those of many other instabilities such as, for example, those observed 
during the growth of binary alloys [H, in flow of water through a rigid hot porous media Q, or in more complex 
geomorphological settings [101 ] , and we structure our model accordingly. 

Katsaros [f| and Woodcock Q attribute the holes from which the stars emanate and the patterns themselves to 
thermal convection patterns within the lake, but do not measure or calculate their nature. However, often the holes 
do not exhibit a characteristic distance between them but rather form from protrusions (e.g. sticks that poke through 
the ice surface) [BJ and stars follow thereby ruling out a convective mechanism as being necessary to explain the 
phenomena. The paucity of literature on this topic provides little more than speculation regarding the puncturing 
mechanism but lake stars are observed in all of these circumstances. Therefore, while hole formation is necessary for 
lake star formation, its origin does not control the mechanism of pattern formation, which is the focus of the present 
work. 



II. THEORY 

The water level in the hole is higher than that in the wet snow-slush-layer Q and hence we treat this warm water 
region as having a constant height above the ice or equivalently a constant pressure head, which drives flow 
of water through the slush layer, which we treat as a Darcy flow of water at 0°C. We model the temperature field 
within the liquid region with an advection-diffusion equation and impose an appropriate (Stefan) condition for energy 
conservation at the water-slush interface. The water is everywhere incompressible. Finally, the model is closed with 
an outer boundary condition at which the pressure head is assumed known. 

Although we lack in-situ pressure measurements, circular water-saturated regions (a few meters in radius) are 
observed around the lake stars. Hence, we assume that the differential pressure head falls to zero somewhere in the 
vicinity of this circular boundary. The actual boundary at which the differential pressure head is zero is not likely 
to be completely uniform (as in Figure 4 of Knight Q) but treating it as uniform is a good approximation in the 
linear regime of our analysis. Finally, we treat the flow as two-dimensional. Thus, although the water in direct 
contact with ice must be at 0°C, we consider the depth-averaged temperature, which is above freezing. Additionally, 
the decreasing pressure head in the radial direction must be accompanied by a corresponding drop in water level. 
Therefore, although the driving force is more accurately described as deriving from an axisymmetric gravity current, 
the front whose stability we assess is controlled by the same essential physical processes that we model herein. Our 
analysis could be extended to account for these three-dimensional effects. 

The system is characterized by the temperature T, a Darcy fluid velocity u, pressure p, and an evolving liquid- 
slush interface a. The liquid properties are k (thermal diffusivity) , Cp (specific heat at constant pressure) and fi 
(dynamic viscosity) and the slush properties are n (permeability), £ (solid fraction) and L (latent heat). We non- 
dimcnsionalizcd the equations of motion by scaling the length, temperature, pressure and velocity with ro, To, po, 
and respectively. Thus, our model consists of the following system of dimensionless equations: 
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f)0 

— + u-V6 = eV 2 8r <r<a(cb,t), (1) 



dt 



with boundary conditions 



a((f>,t) < r < 1, (2) 



p = 1 ri <r < a(<f>,t), (3) 
V 2 p = a((M) <r < 1, (4) 



V • u = n<r < a(<f>,t), (5) 
u| a _=u| a+ r = a(<f),t), (6) 
u = -Wp a((f>,t) < r < 1, (7) 



o = --|vfl r = o(0,t), (8) 
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(9) 

[0 r = 1 

and 

!1 r = n 
1 r = a(0,i), (10) 
r = 1 

where ([lj describes the temperature evolution in the liquid, (j4|) and ([5|) describe mass conservation with a Darcy flow 
(0 in the slush, ([8|) is the Stefan condition, and J9J and (fTU|) are the temperature and pressure boundary conditions, 
respectively (see Figure [2]) ■ Note that ([3]) and ([S]) can both be satisfied since the liquid region has an effectively infinite 
permeability. 

The dimensionless parameters e and S of the system are given by 

e^—, and S = (11) 
uor Op Jo 

which describe an inverse Peclet number and a Stefan number respectively. Because the liquid must be less than 
or equal to 4°C, we make the conservative estimates that T < 4°C, £ > 0.3, and use the fact that L/Cp « 80°C 
from which we see S > 6 > 1. Using k ~ 10 _7 m 2 s _1 , and the field observations of Knight Q to constrain 
(lcm/hr < uq < lOcm/hr) and ro (0.3m < ro < 3m), we find that e < 0.1 -C 1. We therefore employ the quasi- 
stationary (S 3> 1) and large Peclet number (e 1) approximations, and hence equations - (fTO)) are easily 
solved for a purely radial flow with cylindrical symmetry (no <f) dependence) and circular liquid-slush interface. This 
(boundary layer) solution is 

u = ur = — — — --r Ti < r < 1, (12) 
in(a ) r 

ln(r) 

Pb = 7—, — r r > ao, (13) 
ln(ao) 

,. N i(-l/ln(a )+2e) 



% = 1 - ( — ) r < a , (14) 
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FIG. 3: Stability curve: Non-dimensional growth rate a versus non-dimensional wavenumber k! . Scales for the axes are given 
at the upper left (a axis) and the lower right corners (k' axis), a is plotted for the range of plausible ao (dot-dashed blue and 
dashed red curves) and for the approximation (|18p (solid green curve). 



-l/ln(a )+2e ' [ ' 



where equation (|15|) has an approximate implicit solution for ao given by 



^-i §ln(«o) = |- (16) 

We perform a linear stability analysis around this quasi-steady cylindrically symmetrical flow. Proceeding in the 
usual way, we allow for scaled perturbations in 9 and a with scaled wavenumber k' — efc, non-dimensional growth rate 
(T, and amplitudes /(r) and g respectively. Keeping only terms linear in e, 1/S and g, we solve (j4]) subject to (jlOl) . 
substitute into ^ and satisfy © and ([T|). This gives the non-dimensional growth rate (a) as a function of scaled 
wave number (fc'): 



2a \n 2 (a )S V V ~ ' " J \-k'ln(a ) 
Equation (|17p can be approximated in < x < 1 as 

a 



l + 4fc^ln a (ao)-l) I „°°, , -1). (17) 



\n 2 (a )S 



x(l-x), (18) 



where x = — k' ln(ao)/ao- 

The stability curve (|T7|) and the approximation (fT8|) are plotted in Figure [3] The essential features of (|17|) are a 
maximum in the range < k' < ao/\n(ao), zero growth rate at k' — ao/ln(ao) and a linear increase in stability with 
k' for large k'. The long- wavelength cut-off is typical of systems with a Peclet number, here with the added effect of 
latent heat embodied in the Stefan number. This demonstrates the competition between the advection and diffusion 
of heat and momentum (in a harmonic pressure field); the former driving the instability and the latter limiting its 
extent. The maximum growth rate occurs at approximately 



with (non-dimensional) growth rate 



k > « ^ (19) 

max -21n(oo)' 1 ' 



% ■ (20) 

4S , ln 2 (a ) 



FIG. 4: Schematic showing ro, r% , tls and ri. 



Translating (|19p and (|20[) back into dimensional quantities, we find that the most unstable mode has angular size 
given by 

^degrees = ^ ( ~) m ( ~ ) ' ( 21 ) 



uor \a J \a 

and has growth rate given by 



adim ~ 45r l n >o/ao) (ro ' , "' 



III. EXTRACTING INFORMATION FROM FIELD OBSERVATIONS 



Field observations of lake stars cannot be controlled. A reasonable estmate for ro is the radius of the wetted 
(snow) region around the lake stars, and observations @, @, 0] bound the value as 1.5m < ro < 4m. This is simply 
because if there were significant excess pressure at this point then the wetting front would have advanced further. 
However, it is also possible that the effective value of ro, say r$ , is less than this either because the wetted radius 
is smaller earlier in the star formation process or because the ambient pressure level is reached at smaller radii. 
Here, we take ao to be the radius of the roughly circular liquid-filled region at the center of the lake star (re) as 
the best approximation during the initial stages of star formation (see Figure 3]). Field observations show that 
0.1m < re < 0.5m, [1, H, 0] and hence 0.07 < rg/r < 0.15. We note that equations (|21~j) and ((22)) are more sensitive 
to ao/ro than ao or ro independently [13]. With this interpretation of ro we find a reasonable estimate of uo as 
1.4 • 10~ 5 m/s < Uo < 2.8 • 10~ 5 m/s. Using these parameter values, the most unstable mode should have wavelength 
between 8° and 130°. Letting the number of branches be N = 360° /<fideg, then 3 < N < 45 and we clearly encompass 
the observed values for lake stars (4 < N < 15), but note that values (N > 15) are never seen in the field. 

Despite the dearth of field observations, many qualitative features embolden our interpretation. For example, the 
stars with larger values of ao/ro have a larger number of branches. Moreover, for any value of ao/ro, our analysis 
predicts an increase in N with ro and uq. Indeed, uq increases with po (higher water height within the slush layer) 
and n (less well-packed snow). Therefore, we ascribe some of the variability among field observations to variations 
in these quantities (which have not been measured in the field) and the remainder to nonlinear effects. Because the 
dendritic arms are observed long after onset and are far from small perturbations to a radially symmetric pattern, 
as one might see in the initial stages of the Saffman- Taylor instability, the process involves non-linear cooperative 
phenomena. Hence, our model should only approximately agree with observations. Although a rigorous non-linear 
analysis of the long term star evolution process (e.g. 0]) may more closely mirror field observations, the present state 
of the latter does not warrant that level of detail. Instead, we examine the model physics through simple proof of 
concept experimentation described presently. 
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FIG. 5: Typical experimental run where small- scale fingers are present. For scale, the nozzle head has diameter of 5 mm. 




FIG. 6: Typical run where channels form. This picture is taken from the underside. Note: part of the slush broke off when it 
was flipped to image it. The ruler scale is in cm. 



IV. DEMONSTRATING LAKE STARS IN THE LABORATORY 

A 30 cm diameter circular plate is maintained below freezing (« — 0.5°C), and on top of this we place a 0.5 to 1 cm 
deep layer of slush through which we flow 1°C water. Given the technical difficulties associated with its production, 
the grain size, and hence the permeability, of the slush layer, is not a controlled variable. This fact influences our 
results quantitatively. In fourteen runs we varied the initial size of the water- filled central hole (ao), that of the 
circular slush layer (ro), and the flow rate (Q), which determines uq. The flow rate is adjusted manually so that the 
water level (ho) in the central hole remains constant [131 ]. Fingering is observed in every experimental run and hence 
we conclude that fingers are a robust feature of the system. Two distinct types of fingering are observed: small-scale 
fingering (see Figure [5]) that forms early in an experimental run, and larger channel- like fingers (see Figure [S]) that are 
ubiquitous at later times and often extend from the central hole to the outer edge of the slush. Since the channel-like 
fingers provide a direct path for water to flow, effectively shorting Darcy flow within the slush, their subsequent 
dynamics are not directly analogous to those in natural lake stars. However, in all runs, the initial small-scale fingers 
have the characteristics of lake stars and hence we focus upon them. We note that because the larger channel-like 
fingers emerge out of small-scale fingers, they likely represent the non-linear growth of the linear modes of instability, 
a topic left for future study. Finally, we measure the distance between fingers (df), so that for each experiment we 
can calculate uq = Q/(2irroho), 4> ca ic = 't'degrees, from equation (pH]) . and (j) i, s = 180° d f / (irao) , and we can thereby 
compare experiment, theory and field observations. 

In Figure [7] we plot <f> b s versus </> ca / c for the various field observations for which we have estimates of parameters, 
the laboratory experiments described above, and the model [equation (|21[) ]. There is a large amount of scatter in both 
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FIG. 7: Comparison of theory, experiment and field observations. Circles are field observations (cyan = best constrained field 
observation, black = range of plausible field observations), triangles are experimental results (blue upward-pointing triangles 
were unambiguous; red left-pointing triangles have channels but show no clear small-scale fingers, so channel spacing is taken 
for dj; green right-pointing triangles were compromised by the quality of the images). Errors are approximately 0.3 cm, 0.5 
cm, 2 mm, 5 ml/min and 0.2 cm (respectively) for the five measured quantities. All experimental results thus have error bars 
of at least a factor of two in the x-coordinate and 30% in the y-coordinate. Typical error bars are shown on one measurement. 
The solid red line is the theoretical prediction; the dotted green line is the best fit line to the blue triangles. 



the experimental and observational data and the data does not lie on the one-to-one curve predicted by the model. 
However, the experiments are meant to demonstrate the features of the model predictions, and the results have the 
correct qualitative trend (having a best-fit slope of 0.34). We also attempt to find trends in the experimental data 
not represented by the model by comparing y = 4> bs/4>caic vs. various combinations of control parameters (= x) 
including 7*0, clq, ro/ao, r^u®, ro/ao ln(ro/ao) and ln(ro/ao)/(aoUo)- For all plots of y vs. x, our model predicts a zero 
slope (and y-intercept of 1). A non-random dependence of y on x would point to failure of some part of our model. 
Thus, to test the validity of our model, we perform significance tests on all non-flagged data with the null hypothesis 
being a non-zero slope. In all cases, the null hypothesis is accepted (not rejected) at the 95% confidence level. Thus, 
although the agreement is far from perfect, the simple model captures all of the significant trends in the experimental 
data. 



By generalizing and quantifying the heuristic ideas of Knight |5j, we have constructed a theory that is able to 
explain the radiating finger-like patterns on lake ice that we call lake stars. The model yields a prediction for the 
wavelength of the most unstable mode as a function of various physical parameters that agrees with field observations. 
Proof of concept experiments revealed the robustness of the fingering pattern, and to leading order the results also 
agree with the model. There is substantial scatter in the data, and the overall comparison between field observations, 
model and experiment demonstrates the need for a comprehensive measurement program and a fully nonlinear theory 
which will yield better quantitative comparisons. However, the general predictions of our theory capture the leading 
order features of the system. 



V. 



CONCLUSIONS 
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